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Abstract. We study the phase diagram of an SU(3)-symmetric mixture of three- 
component ultracold fermions with attractive interactions in an optical lattice, 
including the additional effect on the mixture of an effective three-body constraint 

induced by three-body losses. Wc address the properties of the system in D > 2 by 
using dynamical mean-field theory and variational Monte Carlo techniques. The phase 
diagram of the model shows a strong interplay between magnetism and superfluidity. 
In the absence of the three-body constraint (no losses), the system undergoes a phase 
transition from a color supcrfiuid phase to a trionic phase, which shows additional 
particle density modulations at half-filling. Away from the particle-hole symmetric 
point the color superfluid phase is always spontaneously magnetized, leading to the 
formation of different color superfluid domains in systems where the total number 
of particles of each species is conserved. This can be seen as the SU(3) symmetric 
realization of a more general tendency to phase-separation in three-component Fermi 
mixtures. The three-body constraint strongly disfavors the trionic phase, stabilizing 
a (fully magnetized) color superfluid also at strong coupling. With increasing 
temperature we observe a transition to a non-magnetized SU (3) Fermi liquid phase. 
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1. Introduction 

Cold atoms in optical lattices provide us with an excellent tool to investigate notoriously 
difficult problems in condensed matter physics [H |2] . Recent progress towards this goal 
is exemplified by the experimental observation of the fermionic Mott insulator [21 S] in a 
binary mixture of repulsively interacting atoms loaded into an optical lattice, and of 
the crossover between Bardeen-Cooper-Schrieffer (BCS) superfiuidity and Bose-Einstein 
condensation (BEC) [HI E], [7j in a mixture of ^Li atoms with attractive interactions. 

At the same time, ultracold quantum gases also allow us to investigate systems 
which have no immediate counterparts in condensed matter. This is the case for 
fermionic mixtures where three internal states u = 1, 2, 3 are used, instead of the usual 
binary mixtures that mimic the electronic spin a =t, i. These multi-species Fermi 
mixtures are already available in the laboratory, where three different magnetic sublevels 
of ^Li [HI [9l [ini [11] or ^^^Yb [12j , as well as a mixture of the two internal states of ^Li with 
a lowest hyperfine state of '^''K[T3] have been successfully trapped. In the case of Alkali 
atoms, magnetic or optical Fano-Feshbach resonances can be used to tune magnitude 
and sign of the interactions in the system, and in the case of Ytterbium or group II 
atoms, it is possible to realise three-component mixtures where the components differ 
only by nuclear spin, and therefore exhibit SU(3) symmetric interactions [HI |15l [16] . 
Moreover, loading these mixtures into an optical lattice would give experimental access 
to intriguing physical scenarios, since they can realize a three-species Hubbard model 
with a high degree of control of the Hamiltonian parameters. 

Multi-species Hubbard models have attracted considerable interest on the 
theoretical side in recent years. First studies were focused on the 5'f/(3)-symmetric 
version of the model with attractive interaction. By using a generalized BCS approach 
[T7] [T8| [T9] , it was shown that the ground state at weak-coupling spontaneously breaks 
the S'f/(3) ® U{1) symmetry down to SU(2) ® f/(l), giving rise to a color superfluid (c- 
SF) phase, where superfluid pairs coexist with unpaired fermions. Within a variational 
Gutzwiller technique [201 EI] the superfluid phase was then found to undergo for 
increasing attraction a phase transition to a Fermi liquid trionic phase, where bound 
states (trions) of the three different species are formed and the S'f/(3)-symmetry is 
restored. More recently [221 ES] , the same scenario has been found by using a self-energy 
functional approach for the half-flUed model on a Bethe lattice in dimension D = oo. It 
was suggested [21] that this transition bears analogies to the transition between quark 
superfluid and baryonic phase in the context of Quantum Chromo Dynamics. 

Both the attractive and the repulsive version of the model was addressed by 
numerical and analytical techniques for the peculiar case of spatial dimension D = 1 
[25] [26] [271 [28] . while Mott physics and instabilities towards (colored) density wave 
formation have been found in the repulsive case in higher dimensions [T71 [29] [30] . 
It is important to mention that substantial differences are expected in the attractive 
case at strong coupling when the lattice is not present [311 132] • Those differences are 
essentially related to the influence of the lattice in the strong coupling limit in the three- 
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body problem, favoring trion formation [331 El] with respect to pair formation in the 
continuum, as was shown in Ref . [321 |35l [36] . 

Here we consider the S'f/(3)-symmetric system in a lattice for D > 2 m the 
presence of attractive two-body interactions by combining dynamical mean-field theory 
(DMFT) and variational Monte Carlo (VMC). We analyze several cases of interest 
for commensurate and incommensurate density. Ground state, spectral, and finite 
temperature properties are addressed. More specifically we focus on the transition 
between color superfiuid and trionic phase and on a better understanding of the 
coexistence of magnetism and superfluidity in the color superfiuid phase already 
predicted in the SU (3) symmetric case [201 [21] but also when the S'f/(3)-symmetry is 
explicitly broken [37] . We show that the existence of a spontaneous magnetization leads 
the system to separate in color superfiuid domains with different realizations of color 
pairing and magnetizations whenever the total number of particles in each hyperfine 
state is conserved. This would represent a special case, due to the underlying SU{^) 
symmetry, of a more general tendency towards phase separation in three-component 
Fermi mixtures. We point out that all this rich and interesting physics arises merely 
from having three components instead of two. Indeed the analogous SU (2) system 
would give rise to the more conventional BCS-BEC crossover, where the superfiuid 
ground state evolves continuously for increasing attraction [3H]- Moreover in the SU{2) 
case superfluidity directly competes with magnetism [39] . 

The case under investigation can be realized with ultracold gases by loading a 
three-species mixture of ^^^Yb [I2] or another group II element such as ^^Sr into an 
optical lattice, or alternatively using ^Li in a large magnetic field. However, some 
reahzations with ultracold atoms are plagued by three-body losses due to three-body 
Efimov resonances [HI [HI [TT], which are not any more Pauli suppressed as in the two- 
species case. The three-body loss properties and their dependence on the magnetic 
field have been already measured for ^Li [SJ [HI [TT], while they are still unknown for 
three-component mixtures of certain group-II elements. Loading a gas into an optical 
lattice could be used to suppress losses, as a large rate of onsite three-body loss can 
prevent coherent tunneling processes from populating any site with three particles [IQ] . 
As proposed in Ref. [^ for bosonic systems, in the strong loss regime a Hamiltonian 
formulation is still possible if one includes an effective hard-core three-body interaction, 
which leads to new interesting physics [41j. The effect of this dynamically generated 
constraint on the fermionic system in D = 1 with attractive interactions was studied in 
Ref. [27], where it was shown that the constraint may help to stabilize the superfiuid 
phase in some regions of the phase diagram. 

For these reasons we also study the effect of including a three-body constraint in the 
model, as representative of an SU{3) symmetric mixture in the strong-loss regime. The 
asymmetric case in the strong loss regime, which is directly relevant for experiments on 
^Li close to a Feshbach resonance, has been already addressed in a separate publication 

m- 

The paper is organized as follows: in the following sections we first introduce 
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the model (Sec. [2]) and then the methods used (Sec. [3]). Later on we present our 
results, focusing first on the unconstrained system (Sec. [i]), for commensurate and 
incommensurate densities and then on the effects of the three-body constraint (Sec. [sj. 
The emergence of domain formation within globally balanced mixtures is discussed in 
detail in Sec. [H Final remarks are drawn in Section [7l 

2. Model 

Three-component Fermi mixtures with attractive two-body interactions loaded into an 
optical lattice are well described by the following Hamiltonian 

'H = - J ^ c\„Cj^„ - ^ ^lani,a + X] X] ^<^^'^i,'^'^i^' + ^ 5Z ^l.^^S.i^-S.i , (1) 
{i,j),a i,(T i a<a' i 

where a = 1, 2, 3 denotes the different components, J is the hopping parameter between 
nearest neighbor sites fJ^a is the chemical potential for the species a and f/g-o-' < 0. 

We introduced the onsite density operators = c\^Cia. The three-body interaction 
term with = oo is introduced to take the effects of three-body losses in the strong loss 
regime into account according to Refs. [271 HO]- V = corresponds to the case when 
three-body losses are negligible. While the model and the methods are developed for 
the general case without SU (3)-symmetry, in this paper we concentrate on the SU (3)- 
symmetric case reflected by species- independent parameters 

f/^^/ = U, /i^ = n. (2) 

In this case the Hamiltonian ([T]) reduces to an SU (3) attractive Hubbard model if = 0. 
Note that the three-body interaction term is a color singlet and thus does not break 
SU (3) for any choice of V. On the basis of previous works, the ground state of the 
unconstrained model is expected to be, at least in the weak coupling regime, a color 
superfiuid, i.e. a phase where the full SU{3) ® U{1) symmetry of the Hamiltonian is 
spontaneously broken to SU{2)^U{1) [l7l[T8]. As shown in [TTlilB], it is always possible 
to find a suitable gauge transformation such that pairing takes place only between two 
of the natural species a, a' and in this paper we choose a gauge in which pairing takes 
place between the species a = 1 and a' = 2 (1 — 2 channel), while the third species 
stays unpaired. Whenever the S'?7(3)-symmetry is explicitly broken, only the pairing 
between the natural species is allowed to comply with Ward-Takahashi identities [37] . 
This reduces the continuum set of equivalent pairing channels of the symmetric model 
to a discrete set of three {mutually exclusive) options for pairing, i.e. 1 — 2, 1 — 3 or 
2 — 3. In this case the natural choice would be that pairing takes place in the channel 
corresponding to the strongest coupling when the mixture is globally balanced. We can 
always relabel the species such that strongest attractive channel is the channel 1 — 2. 
Other pairing channels can be studied via index permutations of the species. Therefore 
the formalism developed here is fully general and includes both the symmetric and non- 
symmetric case, while only in the S'f/(3)-symmetric case our approach corresponds to a 
specific choice of the gauge. 
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3. Methods 

In order to investigate the model in Eq. ([T]) in spatial dimensions D > 2 we use a 
combination of numerical techniques which have proven to give very consistent results 
for the non-symmetric case |12]. In particular, we use dynamical mean- field theory 
(DMFT) for L> > 3 and variational Monte Carlo (VMC) for D = 2. DMFT provides 
us with the exact solution in infinite dimension and a powerful (and non-perturbative) 
approach in D = 3, which has the advantage of being directly implemented in the 
thermodynamic limit (without finite size effects). VMC allows us to incorporate also 
the effect of spatial fluctuations which are not included within DMFT, even though the 
exponential growth of the Hilbert space limits the system sizes that are accessible. 

3.1. DMFT 

Dynamical mean-field theory (DMFT) is a non-perturbative technique based on the 
original idea of Metzner and VoUhardt who studied the limit of infinite dimension of 
the Hubbard model In this limit, the self-energy Zl(k, cu) becomes momentum 

independent S(k, w) = S(a;), while fully retaining its frequency dependence. Therefore 
the many-body problem simplifies significantly, without becoming trivial, and can be 
solved exactly. In this sense DMFT is a quantum version of the static mean-field theory 
for classical systems, since it becomes exact in the same limiting case [D = oo) and 
can provide useful information also outside of this limit, fully including local quantum 
fluctuations. In 3D, assuming a momentum independent self-energy, has proved to be 
a very accurate approximation for many problems where the momentum dependence 
is not crucial to describe the physics of the system such as the Mott metal-insulator 
transition [H] where the frequency dependence is more relevant than the k dependence. 

3.1.1. Theoretical setup for SU{3) model with spontaneous symmetry breaking - In 
this work, we generalize the DMFT approach to multi-species Fermi mixtures in order 
to describe color superfiuid and trionic phases, which are the expected phases occurring 
in the system. The theory can be formulated in terms of a set of self-consistency 
equations for the components of the local single-particle Green function G on the lattice. 
Since we are dealing here with superfiuid phases involving also anomalous components 
of the Green function, we use a compact notation in terms of mixed Nambu spinors 
ip = (ci,C2,C3), where we already assumed that pairing takes place only between the 
first two species, as explained in the previous section, and we omit the subscript i 
(spatially homogeneous solution). We reiterate that this specific choice is valid without 
loss of generality in the SU (3)-symmetric model, and has the same status as fixing the 
phase of a complex condensate order parameter in theories with global phase symmetry. 



Magnetism and domain formation in SU (3) -symmetric multi-species Fermi mixtures 6 
The local Green function G{iUn) in Matsubara space then has the form 

/ GiilUJn) F{lUJn) \ 

G{lUJn)= F*{-tUn) -GlitUn) , (3) 

\ G^itUn) J 

where G^ir) = —{TrCa{T)c\) and F{t) = —{TrCi{T)c2) are respectively the normal 
and anomalous Green functions in imaginary time, and G(j{iujn) = dTG(j{T)e'^^^'^ 
and F{iuJn) = /g'^ F(r)e*^"'^ are their Fourier transforms in Matsubara space, where 
Un = (2n + l)7rT (ks = 1). 

In practice the original lattice model ([T| in the DMFT approach can be mapped, 
by introducing auxiliary fermionic degrees of freedom a|^, a^^, on a Single Impurity 
Anderson Model (SIAM), whose Hamiltonian reads 

T-LsiAM = Ugg'Tign^' + Vnin2n3 - y^^figU^ (4) 

ct<ct' (T 

+ [eiattlai^ + Vi^ [clai^ + /i.e.)] +YWi al^aj^ + h.c. , 

Icr I 

where the Anderson parameters eia, Via, Wi have to be determined self-consistently. Self- 
consistency ensures that the impurity Green function of the SIAM is identical to the local 
component of the lattice Green function. The components of the non-interacting Green 
function for the impurity site, which represent the dynamical analog of the Weiss field 
in classical statistical mechanics, can be expressed in terms of the Anderson parameters 
as 

71^ y,iQ,2 + 



1T-S T/2 A* 

1=1 ^''3 

where = —ii^n + ^la- The self-consistency equations for the local Green functions 
now have the form 

G{icOn) = j^Yl '^«) = / deD{e)&^'\e, zun) , (9) 

k ^ 

where M is the number of lattice sites, (^'"^^(k, zu;„,) = (^'""(ek, ^i^n) is the lattice 
Green function within DMFT and D[e) is the density of states of the lattice under 
consideration. The independent components of G'"**(k, icjn) have the form 

Qlatt C2 ~ Qg\ 
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Glatt 
2 



-tlatt 



(C2 - £k)(Cl* - ^k) + ^sc{-i^n 



)j:*g^{iuJr^ 



sc 



^latt 



(Ci-^k)(C2*-^k) + s|^( 
1 



where 



Cs - £k ' 

iUn + ficT -^c 



(11) 

(12) 
(13) 



and the self-energy can be obtained by the following 



local Dyson equation J]{iUn) = QatuI^ 

( Si(iw„) 



lUJr. 



T,{iU]n) 



Ssc(«tJn) 



where 
\ 



V 



(14) 



MlCOn) ) 

Once a self-consistent solution has been obtained, the impurity site of the SIAM 
represents a generic site of the lattice model under investigation. Therefore several 
static thermodynamic quantities can be directly evaluated as quantum averages of the 
impurity site. As evident from the previous equations, DMFT is explicitly formulated 
in a grand canonical approach where the chemical potentials [i^ are given as input and 
the onsite densities n^, = (cJ.Co-) are calculated. 

3.1.2. Calculated observables and numerical implementation - To characterize the 
different phases, we evaluated several static observables such as the superfluid (SF) 
order parameter P = (C1C2), the average double occupancy d^^-' = ('^o-^o-') and the 
average triple occupancy t = {nin2nz). As suggested in Refs. [TTl [181 Elj, in order 
to gain condensation energy in the c-SF phase, it is energetically favorable to induce 
a finite density imbalance between the paired species (1 — 2 in our gauge) and the 
unpaired fermions. To quantitatively characterize this feature we introduce the local 
magnetization 



m 



ni2 — where ni2 = ni = 712 . 



(15) 



From the normal components of the lattice Green functions in Eqs. (10), and 



(13) we can extract the DMFT momentum distribution 
n,(k)=Tj]G'|,»**(k,zu;„)e- 



and the average of kinetic energy per lattice site 

K = j^^i^knaiX) = ^ I de D{e) e n^{e). 



(16) 



(17) 



It is evident from the expression of (j'"**(k, iw^) given in Eqs. (10), (11) and (13) that 
no^(k) only depends on the momentum k through the free-particle dispersion of the 
lattice at hand. The internal energy per lattice site E can then be obtained as 

E = K ^ Vpot, where Vpot = ^ ^^rf^a' (18) 
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is the average potential energy per lattice site. 

Solving the DMFT equations is equivalent to solving a SIAM in presence of a bath 
determined self-consistently. We use Exact Diagonalization (ED) [15], which amounts 
to truncating the number of auxiliary degrees of freedom a|^, a^^ in the Anderson model 
to a finite (and small) number Ng — 1. In this way the size of the Hilbert space of the 
SIAM is manageable and we can exactly solve the Anderson model numerically. Here 
we would like to point out that this truncation does not reflect the size of the physical 
lattice but only the number of independent parameters used in the description of the 
local dynamics. Therefore we always describe the system in the thermodynamic limit 
(no finite-size effects). We use the Lanczos algorithm [IB] to study the ground state 
properties (up to Ng = 7) and full ED for finite temperature (up to A^s = 5). Due to 
the increasing size of the Hilbert space (cr = 1,2,3 instead of a =t;i) in the multi- 
component case the typical values of Ng which can be handled sensibly is smaller than 
the corresponding values for the SU{2) superfiuid case. However, in thermodynamic 
quantities, we found indeed only a very weak dependence on the value of Ng and the 
results within full ED at the lowest temperatures are in close agreement with T = 
calculations within Lanczos. 

A definite advantage of ED is that it allows us to directly calculate dynamical 
observables for real frequencies without need of analytical continuation from imaginary 
time. In particular, we can directly extract the local single-particle Green function 
G'o-(w) and the single-particle spectral function 

p^ico) = --ImG^ico + iO+). (19) 

71 

3.2. Variational Monte-Carlo 

The variational Monte Carlo (VMC) techniques described in this subsection can be 
used to calculate the energies and correlation functions of the homogeneous phases at 
T = in a canonical framework. The basic ingredients of the VMC formalism are the 
Hamiltonian and trial wavefunctions with an appropriate symmetry. In principle, the 
formalism presented here can be applied to any dimension, even though here we use it 
specifically to address the system on a two-dimensional square lattice. 

The canonical version of Hamiltonian ([T| for three-components fermions with 
generic attractive interactions is given by 

n = -jJ2 ^34aC„.^3 + J2 Uaa'n,,,ni^„>, (20) 
{i,j),(J i,a<cr' 

where the three-body constraint is imposed by using the projector V3 = Ylii^ ^ 
^i,i^j,2'^i,3) and in the unconstrained case we set V3 equal to the identity. 

Practical limitations do not permit a general trial wave function equally accurate 
both for the weak- and the strong-coupling limit. Due to this reason we introduce 
different trial wavefunctions for different coupling regimes. 

In the weakly interacting limit, which we operatively define as \U(j„r \ < 4 J = W/2, 



we use the full Hamiltonian ( 20 ) along with the weak-coupling trial wavefunction defined 
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in the next subsection. Here W = 2DJ is the bandwidth. At strong-couphng this 
wavefunction results in a poor description of the system. In order to gain insight into 
the strong couphng regime, we derive below a perturbative Hamiltonian to the second 
order in J/Uaa', which we will combine with a strong-coupling trial wavefunction. Again 



the strong-coupling wavefunctions are incompatible with the Hamiltonian (20), as will 
be clarified below. We can therefore address confidently both limits of the model while 
at intermediate coupling we expect our VMC results to be less accurate. 

3.2.1. Strong coupling Hamiltonian, constrained case - In order to derive a 
perturbative strong-coupling Hamiltonian for the constrained case we make use of the 
Wolff-Schrieffer transformation [TT] 

npert = Voe'^ne-'^VD (21) 

and keep terms up to the second order in J/Uaa'- In the expression above, Vd is the 
projection operator to the Hilbert subspace with fixed numbers of double occupancies in 



each channel {N^', Nj^,N^^), and e*"^ is a unitary transformation defined in Appendix 
[a} So we obtain the perturbative Hamiltonian (see [Appendix A ), which reads: 



{i,j)cr {j',i);{i,j);cr<(7' 

-J' E 77^4.v//,.'4rf.,..' + V + 0(-^) , (22) 

(i,i'>;(ij>;-<-' 

where V = XI Uaa'^i^a^iy. Here we define double occupancy operators as dl^^, = 

cl a^iyhiyi {hi^u = 1 — Ui^a) and single occupancy operators as //^ = hiyhiyc] ^ . 

For the case where the S'f/(3)-symmetry is restored (f/o-o-' = U), the perturbative 
Hamiltonian can be written in a compact notation 

n,ert = v-jY1 [fU^''^ + 4,A,<] -jf H 4,.A<^flA,<^ (23) 

where the double occupancy operator is now defined as rfj^. = cJo-(^«,o-''"'j,o-" + ^j,o-"'^j,cr')- 
Now, rather than conserving the number of double occupancies N^'^' in each 
channel, only the total number A^'^^o = + N^^ + Nf" is conserved due to the 



S'f/(3)-symmetry. Indeed Eq. (23) contains terms where the tightly bound dimers 
are allowed to change the composition through second order processes. Thus, the 
SU (3)-symmetric case, in contrast to the case with strongly anisotropic interactions 
is qualitatively different from the Bose-Fermi mixture, because the bosons - tightly 
bound dimers - can change composition as described above, while such a process was 
not allowed in the case of the strong anisotropic interactions. We also notice that the 
last of the ~ J^/t/ terms contributes only when N^fl < N/2. 
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3.2.2. Strong coupling Hamiltonian, unconstrained case - Without the 3-body 
constraint three fermions with different hyperfine states can occupy the same lattice 
site and we expect them to form trionic bound states at sufficiently strong coupling. 
Correspondingly the many-body system should be in a trionic phase with heavy trionic 
quasiparticles, as mentioned in previous studies [TTl HSl 123] • Therefore we expect 
that our perturbative approach can provide a description of the trions in the strong 
coupling limit. 

First we consider the extreme case J = 0. In this limit formation of local trions 
takes place, i.e each site is either empty or occupied by three fermions with different 
hyperfine spins. Their spatial distribution is random, because any distribution of trions 
will have the same energy. For finite J with J <^ |f/cr.(T'| the hopping term can break a 
local trion, but this would result in a large energy penalty. 

According to perturbation theory up to third order we could have two different 
contributions: (i) one of the fermions hops to one of the neighboring sites and returns 
back to the original site (second order perturbation), (ii) all three fermions hop to the 
same nearest neighbor site (third order perturbation). As we show below, due to the 
first process there is an effective interaction between trions on nearest neighbor sites. 
Also due to this process the onsite energy has to be renormalized. The second process 
(ii) describes the hopping of a local trions to a neighboring site. 

After straightfroward calculations (see Appendix A) we obtain that the effective 
interaction between two trions on neighboring sites is 

l4„^AB.-AE„^-(j^ + ^^ + j^). (24, 

For the SU (3)-symmetric case this expression is simplified and we obtain 

Therefore the nearest-neighbour interaction between trions is repulsive in the SU{?))- 
symmetric case. 

For the hopping coefficient we obtain 

^ ^ (f/... + t/<,..)(f/<.." + f/.v") ^^^^ 

cr,(7 

where a, o' and o" are different from each other in the sum. 

In the SU (3)-symmetric case, the expression again simplifies to 
3 73 

= ^ (27) 



So we obtain the following effective Hamiltonian 

^e// = -Jeff tit, + V^ff ^ njnj . (28) 

Here t] is the creation operator of a local trion at lattice site i and nj = t\ti is the 
trionic number operator. Because the effective hopping of trions results from a third 
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order process and the interaction from second order, more precisely Jeff = J/\U\ ■ Veff, 
the effective trion theory is interaction dominated. Since the interaction describes 
nearest-neighbour repulsion, the strong coupling limit clearly favors a checkerboard 
charge density wave ground state at half-filling [H which we will discuss in more detail 
in Sec. HI 



3.2.3. Trial wavefunctions: In order to describe a normal Fermi liquid phase without 
superfluid pairing, we use the following trial wavefunction 

\NFF) = JVsPdII n i^lO)' (29) 

where |0) is the vacuum state and ^k.o- = —2J{cos{kx) -\-cos{ky)) for a 2D square lattice 
with only nearest-neighbor hopping. The dependence on the densities is included in 
the value of the non-interacting Fermi energy eF,a- The wavefunction above has no 
variational parameters except for the choice of Jastrow factor 

^ f exp(z/3 ni,inj^2^j,3) unconstrained case, weak coupling 

exp(i/c Z](ij)('^i,i^*,2^i,3 + ni^ini,3nj^2 + ni,2^i,3'^i,i)) constrained case 

which takes into account the effect of the interaction. Here and are variational 
parameters and is summation with nearest neigbours. The weak-coupling version 

of the wavefunctions presented in this part is obtained by setting Vd equal to unity. 

We also consider the broken symmetry SU{2) (g) U (1) phase with s-wave pairing in 
the 1 — 2 channel, whose trial wavefunction is given by 

|c - SF) = JV.Vd n [«k + v^.c\A,^ n ^i'.slO) ' (31) 

where = ^ ^1 + (ck — i^)/ \/ (^k — A^)^ + (^''(k))^ j and v\ = \ — m^. In this case, 
in addition to the Jastrow factor JT", we have fi and Aq as additional variational 
parameters. The s-wave gap function A*(k) = Aq has no k dependence. This 
parametrization of A'*(k) leads upon Fourier transform to a singlet symmetric pairing 
orbital 0'^(ri,r2) = 0'^(r2,ri). 

In practice the optimization parameter Aq depends on the density n as well as on 
the coupling strength \J . Also, even at the same coupling strength f/, the Aq can be 
qualitatively different for the weak and the strong coupling ansatz (in the intermediate 
regime f/ ^ —5 J). On the other hand, the parameter fx depends mostly on n (and only 
weakly on t/). The general tendency we observe is that Aq is suppressed beyond the 
filling density n > 1 in the presence of the constraint. Within a BCS mean-field theory 
approach, the condensation energy Econd is easily related to the order parameter Aq, 
being Ecand oc Aq. We however calculate it explicitly from the definition by comparing 



I Despite our fermions are charge neutral, we use sometimes the expression charge density wave in 
analogy with the terminology commonly used in condensed matter physics. 
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the ground state energies of the normal and the superfluid phases for the same density. 
Therefore we define 

Econd = EnfF — Ec-SF , (32) 

where 

Ej^FF = {NFF\n\NFF)/{NFF\NFF) , (33) 

E^_SF = (c - SF\n\c - SF)/{c - SF\c - SF) . (34) 

We also calculate the order parameter P that characterizes the superfluid correlation 
by considering the long range behavior of the pair correlation function 



^".l-^W- (35) 



M 

3 



where B\ = cl^icl^2 M is the total number of the lattice sites. 

Finally, in order to describe the trionic Fermi liquid phase we can use the following 
trial wavefunction 

\Trzon) = Jt JJ 4|0), (36) 

ek<eF 

In this case the Jastrow factor 

Ji = ^exp(-z/i J]nfnJ). (37) 

<«j> 

Here z/f is a variational parameter and summation over nearest neigbours. 



4. Results: SU{3) attractive Hubbard Model 

We first consider the SU{3) attractive Hubbard model described by the Hamiltonian 
([T]) with ^ = 0. In a physical realization with ultracold gases in optical lattices, this 
corresponds to a situation where three-body losses are negligible. In order to address 
the effects of dimensionality and of particle-hole symmetry, we analyze several cases of 
interest, namely (i) an infinite-dimensional Bethe lattice in the commensurate case (half- 
filling), (ii) a three-dimensional cubic lattice and (iii) a two-dimensional square lattice, 
the latter two in the incommensurate case. In order to simplify comparison of results 
on different dimensions, we rescaled everywhere the energies by the bandwidth W of 
the specific lattice under consideration. For a Bethe lattice in D = oo the bandwidth 
is related to the hopping parameter hj W = 4J, while for a D dimensional hypercubic 
lattice it is = ADJ. 



Bethe lattice at half-filling 

We first consider the infinite dimensional case, for which the DMFT approach provides 
the exact solution of the many-body problem whenever the symmetry breaking pattern 
of the system can be correctly anticipated. For technical reasons we consider here the 
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Bethe lattice m D = oo, which has a well defined semicircular density of states, given 
by the following expression 



Die) = ^vW2F^ (38) 

The simple form of the self-consistency relation for DMFT on the Bethe lattice 
introduces technical advantages, as explained below. Moreover, we can directly compare 
our results with recent calculations for the same system within a Self-energy Functional 
Approach (SPA) [221 [23]. 

In the absence of three-body repulsion, the Hamiltonian ([T]) is particle hole- 
symmetric whenever we choose ji = U . In this case the system is half-filled, i.e. no- = | 
for all of (J and n = = 1.5. 

We first consider the ground state properties of the system which we characterize 
via the static and dynamic observables defined in Sec. |3} For small values of the 
interaction (1^71 <^ W), we found the system to be in a color-SF phase, i.e. a phase 
where superfiuid pairs coexist with unpaired fermions (species 1-2 and 3 respectively in 
our gauge) and the superfiuid order parameter P (plotted in Fig. [l] using green triangles) 
is finite. This result is in agreement with previous mean field studies [171 HE], as expected 
since DMFT includes the (static) mean-field approach as a special limit, and with more 
recent SFA results [221 123] . By increasing the interaction \U\ in the c-SF phase, P 
first increases continuously from a BCS-type exponential behavior at weak-coupling to 
a non-BCS regime at intermediate coupling where it shows a maximum and then starts 
decreasing for larger values of \U\. This non-monotonic behavior is beyond reach of 
a static mean- field approach and agrees perfectly with the SFA result [221 123]. As 
explained in the introduction, the spontaneous symmetry breaking in the c-SF phase is 
generally expected [201 ED EH HH] to induce a population imbalance between the paired 
channel and the unpaired fermions, i.e. a finite value of the magnetization m in Eq. 
(15). It is however worth pointing out that, due to particle-hole symmetry, the c-SF 
phase at half-filling does not show any induced population imbalance, i.e. m = for all 
values of the interaction strength. As discussed in the next subsection, the population 
imbalance is indeed triggered by the condensation energy gain in the paired channel. 
This energy gain, however, cannot be realized at half-filling where the condensation 
energy is already maximal for a given U. 

Further increasing \U\, we found P to suddenly drop to zero at \U\ = Uc,2 ~ 0.45iy, 
signaling a first order transition to a non-superfiuid phase. This result is in good 
quantitative agreement with the SFA result in Ref. [221 [23], where a first order transition 
to a trionic phase was found, while a previous variational calculation found a second 
order transition [201 EI] • 

In this new phase we were not able to stabilize a homogeneous solution of the 
DMFT equations with the ED algorithm [15] . Such a spatially homogeneus phase would 
correspond to having identical solutions, within the required tolerance, at iteration n 
and n -|- 1 of the DMFT self-consistency loop. In the normal phase (|t/| > Uc,2) instead, 
we found a staggered pattern in the solutions and convergence is achieved if one applies 
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Figure 1. (Color online) SF order parameter P (green triangles) and CDW order 
parameter C (blue squares) plotted as a function of interaction strength on 
the Bethe lattice in the limit _D — ?> cxo at half-filling and T = 0. C> (empty squares) 
and C< (full squares) correspond to calculations starting from a superfluid or a trionic 
charge density (t-CDW) initial conditions, and similar for P> (note that P< is always 
vanishing and therefore not shown). In the inset we compare the ground state energies 
of the c-SF and t-CDW phases. (Unconstrained, i.e. V = Q) 



a staggered criterion of convergence by comparing the solutions in iteration n and 
n + 2. This behavior is clearly signalling that the transition to a non-superfluid phase is 
accompanied by a spontaneous symmetry breaking of the lattice translational symmetry 
into two inequivalent sublattices A and B. In a generic lattice a proper description of this 
phase would require solving two coupled impurity problems, i.e. one for each sublattice, 
and generalizing the DMFT equations introduced in the previous section. In the Bethe 
lattice instead the two procedures are equivalent [§| 

In the new phase the full S't/(3)-symmetry of the hamiltonian is restored and 
we identify it with as a trionic Charge Density Wave (t-CDW) phase. In order to 
characterize this phase, we introduce a new order parameter which measures the density 
imbalance with respect to the sublattices A (majority) and B (minority), i.e. 

C = ^\nA - ubI (39) 

where ua = n^j^A and rib = n^j^B for all a and C = in the c-SF phase because the 
translational invariance is preserved. The evolution of the CDW order parameter C in 
the t-CDW is shown in Fig. [T] using blue squares. At the phase transition from c-SF 
to t-CDW phase, P goes to zero and C jumps from zero to a finite value. Then C 
increases further with increasing attraction \U\ and eventually saturates at C = 1/2 
for \U\ — ?► oo. Motivated by these findings, we considered more carefully the region 
around the transition point. Surprisingly we found that upon decreasing \U\ from 
strong- to weak-coupling the t-CDW phase survives far below Uc,2 down to a lower 



§ On the Bethe lattice the sublattices A and B are completely decoupled from each other at a given 
step n. 
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Figure 2. (Color online) Sketch of the spatial arrangement of trions in the trionic 
CDW phase. 

critical value Ud — 0, revealing the existence of a coexistence region in analogy with 
the hysteretic behavior found at the Mott transition in the single band Hubbard model 
|44j . In the present case, however, we did not find any simple argument to understand 
which phase is stable and had to directly compare the ground state energy of the two 
phases in the coexistence region to find the actual transition point. In the Bethe 
lattice, the kinetic energy per lattice site K in the c-SF and t-CDW phases can be 
expressed directly in terms of the components of the local Green function G{iUn), which 
is straightforwardly determined by DMFT. The potential energy per lattice site V is 
given by Vt-cnw = ^^^^^^t^^ where the index indicates the sublattice. By generalizing 
analogous expressions valid in the SU{2) case [SHlEniEI], we obtain 

K^.SF = TY,{W/Af[Y^Gl{lUn) - F\lUJn)] (40) 
n a 

and 

^t-CDW = 

Results shown in the inset of Fig. [T] indicate that the t-CDW phase is stable in a 
large part of the coexistence region and that the actual phase transition takes place 
at \U\ = Uc ^ 0.2W. The good agreement between our findings and the SFA results 
in Ref. |22l ES] concerning the maximum value of the attraction Uc2 where a c-SF 
phase solution is found within DMFT would suggest that this value is indeed a critical 
threshold for the existence of a c-SF phase. On the other hand we also proved that the 
c-SF phase close to Uc2 is metastable with respect to the t-CDW phase and therefore the 
existence of the threshold could equally results from an inability of our DMFT solver to 
further follow the metastable c-SF phase at strong coupling. The disagreement between 
our findings and Ref. [221 123] for what concerns the existence of CDW modulations in 
the trionic phase is clearly due to the constraint of homogeneity imposed in the SFA 
approach of Ref. [221 123] in order to stabilize a (metastable) trionic Fermi liquid instead 
of the t-CDW solution. In our case, this was not an issue due to the fact that the 
iterative procedure of solution immediately reflects the spontaneous symmetry breaking 
of the translational invariance and does not allow for the stabilization of an (unphysical) 
homogeneous trionic Fermi liquid at half-filling. 

On the other hand, the necessary presence of CDW modulation in the trionic phase 
at half-filling, at least in the strong-coupling limit, can be easily understood based 
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on general perturbative arguments. Indeed, as pointed out in Sec. [3} in the strong- 
coupling trionic phase where J/\U\ ^ 1, the system can be described in terms of an 
effective trionic Hamiltonian (28). In this Hamiltonian the effective hopping Jg// of 
the trions is much smaller than the next-neighbor repulsion Veff between the trions 
Jeff = ^ Veff = Ij^Tj-- Due to the scaling of the hopping parameter required to 
obtain a meaningful limit D — )■ oo, i.e. J — )■ J j \fz where z is the lattice connectivity, 
one finds Je// — t- in this limit, i.e. the trions become immobile while their next- 
neighbor interaction term survives. In this limit, the Hamiltonian is equivalent to an 
antiferromagnetic Ising model (spin up corresponds to a trion and spin down corresponds 
to a trionic-hole). At half-filling, clearly the most energetically favorable configuration 
is therefore to arrange the trions in a staggered configuration |52]. Moreover, due to 
quantum fluctuations, if we decrease the interaction starting from very large |?7|, the 
spread of a single trion (which is proportional to J'^ /U) increases and it is not a local 
object any more. In this case the trionic wave-function extends also to the nearest 
neighboring sites [3l], as sketched in Fig. [2| This interpretation is in agreement with 
the observed behavior of the CDW order parameter C in Fig. [Tj Indeed, at large 
|f/|, C asymptotically rises to the value C = 1/2, corresponding to the fully local 
trions in a staggered CDW configuration. The presence of the CDW also explains the 
anomalously large value of residual entropy per site Sres = ln2 found when imposing 
a homogeneous trionic phase as in Ref. [22l[23]. At strong-coupling in finite dimensions, 
even though the trions have a finite effective hopping Jeff, one would still expect that the 
augmented symmetry at half-filling favors CDW modulations with respect to a trionic 
Fermi liquid phase. In D = 1, 2 it is indeed known [TTl [27] that the CDW is actually 
stable with respect to the SF phase at half-filling for any value of the interaction, in 
contrast to the SU{2) case where they are degenerate [3H]. Our results prove that in 
higher spatial dimensions this is not the case and there is a finite range of attraction at 
weak-coupling, where the c-SF phase is actually stable. 

Further confirmation of the physical scenario depicted above is provided by the 
analysis of the single-particle spectral function p^. in the c-SF and t-CDW phases shown 
in Fig. |3} In the c-SF phase (Fig. [3]^a)), the spectrum shows a gapless branch due to 
the presence of the third species which is not involved in the pairing, while the spectral 
function for species 1 (2 is identical) shows a gap. The situation is totally different in the 
t-CDW phase (Fig. [sj^b)), where the spectral functions for the three species are identical 
but the lattice symmetry is broken into two sublattices. If we plot the spectral functions 
for the two sublattices (corresponding to two successive iterations in our DMFT loop) 
a CDW gap is visible. We would like to note that the sharply peaked structure of the 
spectrum is due to the finite number of orbitals in the ED algorithm. However, the 
size of the gap should not be affected significantly by the finite number of orbitals. 
Interestingly for \U\ = 0.751^ the size of the energy gap Agap ~ is in very close 
agreement with the value obtained within SFA for the same value of the interaction 
j , indicating that the gap most likely is only weakly affected by CDW ordering. 
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Figure 3. (Color online) Single particle spectral function for the Bethe lattice with 
£) ^ oo at half-fining and T = for (a) the c-SF phase at \U\/W = 0.35 and (b) the 
t-CDW phase at |f/|/VF = 0.75. In the subfigure (a) we plotted pi(a;) (red/dashed 
line) together with —p-i{uj) (green/solid line) to emphasize the different behavior in the 
paired channel and for the unpaired species. The inset shows the low-energy region 
and the c-SF gap. The subfigure (b) shows the spectral function for sublattices A 
(red/dashed line) and B (green/solid line) and the gap in the trionic CDW phase. 
(Unconstrained, i.e. V = Q) 
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Figure 4. (Color online) (a) c-SF order parameter P and (b) CDW amplitude C 
as a function of temperature T/W on the Bethe lattice with D — ^ oo at half-filling. 
Different lines correspond to different values of the interaction. (Unconstrained, i.e. 
F = 0) 



In order to characterize the system at finite temperature, we studied the evolution 
of the SF order parameter P as a function of temperature in the c-SF phase for different 
values of the coupling (Fig. |4]^a)) and analogously for the CDW order parameter C in 
the t-CDW phase (Fig. |4|^b)). The superfluid-to- normal phase transition at T^^{U) is 
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Figure 5. (Color online) Single particle spectral function on the Bethe lattice with 
-D — > oo at half-filling for |?7|/PF = 0.375. Different colors correspond to different 
values of temperature. (Unconstrained, i.e. V — 0) 



also mirrored in the behavior of the spectral function for increasing temperature. The 
results shown in Fig. [5] indicate that the superfluid gap in the spectral function closes 
for T > T^^{U), signaling the transition to a normal homogeneus phase without CDW 
modulations. 

At finite temperatures we also found a coexistence region of the trionic CDW wave 
phase and the color superfluid or normal homogeneous phases in a finite range of the 
interaction U {Ud < \U\ < Uc2 at T = 0). We however leave a thorough investigation of 
the stability range of the t-CDW phase at finite temperature to future study, together 
with its dependence on the distance from the particle-hole symmetric point and on the 
dimensionality. Due to this coexistence region, we define the two critical temperatures 
T^^{U) and T^^^{U) plotted in the phase diagram in Fig. |6| where P{T)\u and 
C{T)\u vanish respectively above the c-SF phase and t-CDW phase. In agreement with 
the results obtained within SFA [221 123], we also found that the critical temperature 
T^^{U) has a maximum at T^^ /W ^ 0.025 for \U\/W = 0.4. This is also in qualitative 
agreement with the SU{2) case [SSj, where the critical temperature has a maximum at 
intermediate couplings. Due to the presence of the CDW modulations in the trionic 
phase which are ignored in Ref . |221 123] , we found also a second critical temperature 
j^CDW ^]2ere charge density wave modulations in the trionic phase disappear. 

4-2. Incommensurate density 

In this section we consider the system for densities far from the particle-hole 
symmetric point. Specifically we investigate, using VMC and DMFT respectively, the 
implementation of the model ([T]) on a simple-square (cubic) lattice in 2D (3D) with 
tight-binding dispersion, i.e. ek = —'2JJ2i=xy{ z) cos{kia), where a is the lattice spacing. 
In particular, we will find that away from the particle-hole symmetric point in the c-SF 
phase, the superfluidity always triggers a density imbalance, i.e. a magnetization. 
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Figure 6. (Color online) Phase diagram of the unconstrained model {V = 0) on 
the Bethe lattice with I? — > oo at half-filling as a function of the temperature T and 
interaction strength \U\. The blue solid line T^^ marks the transition between c-SF to 
a normal phase, while the orange dashed line f^^^ marks the disappearance of CDW 
modulations in the trionic phase. The dashed vertical lines mark the boundaries of 
the coexistence region between the c-SF phase and the t-CDW phase at T = 0. 
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In order to address this feature quantitatively, we studied the system by adjusting 
the chemical potential /i in order to fix the total density n = n^, allowing the system 
to adjust spontaneously the densities in each channel. Due to the spontaneous symmetry 
breaking of the SU{3) symmetry of the Hamiltonian in the color superfluid phase, it 
is indeed possible that, for a given chemical potential /^i = /i2 = yUs = the particle 
densities for different species may differ. If such a situation occurs, the systems shows 
a finite onsite magnetization m. As a more technical remark, we add that the choice of 
pairing channel, as explained in Sec. 3.1.1 , is done without loss of generality: A specific 
choice will therefore determine in which channel a potential magnetization takes place, 
but not infiuence its overall occurrence. Here, since we fix the pairing to occur between 
species 1 and 2, we found a nonzero value of the magnetization parameter m = — n^, 
where = rii = n2. Therefore the paired channel turns out (spontaneously) to be 
fully balanced, while there is in general a finite density imbalance between particles in 
the paired channel with respect to the unpaired fermions. 

The implications of the results presented in this subsection and in Sec. |5] for cold 
atom experiments, where the total number of particles of each species N„ = rij o- is 
fixed, will be discussed in Sec. |6| Combining the grand canonical DMFT results with 
energetic arguments based on canonical VMC calculations, we show that the system is 
generally unstable towards domain formation. 

We first consider in Fig. [7] how the ground state properties of the 3D system evolve 
by fixing the coupling at |f/|/iy = 0.3125, where the system is always found to be in the 
c-SF phase for any density. We consider only densities ranging from n = to half-filling 
n = 1.5. The results above half-filling can be easily obtained exploiting a particle-hole 
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Figure 7. (Color online) c-SF order parameter P (green circles), magnetization m 
(red squares) and average triple occupancy t = (nin2n3) (violet diamonds) plotted as 
a function of the total density n per lattice site for |C//iy| = 0.3125 and T = on the 
cubic lattice in D — 3. The inset shows the behavior of the magnetization in detail. 
(Unconstrained, i.e. y = 0) 

transformation. In particular one easily obtains 

P{n) = P{3 — n) and m{n) = —m{3 — n) , (42) 
t{n) = -t{3-n)+n-2 + d{3-n), (43) 

where t and d are the average triple and double occupancies. The superfluid order 
parameter P increases (decreases) with the density for n < 1.5 {n > 1.5) and is maximal 
at half-filling. The average triple occupancy is instead a monotonic function of the 
density. Below half-filling, the magnetization m first grows with increasing density, then 
reaches a maximum and eventually decreases and vanishes at half-filling in agreement 
with the findings in the previous subsection. This means that in the c-SF phase for a 
fixed value of the chemical potential /x the system favors putting more particles into the 
paired channel than into the unpaired component. For n > 1.5 the effect is the opposite 
and m < 0. This behavior can be understood by considering that the equilibrium value 
of the magnetization results from a competition between the condensation energy gain in 
the paired channel on one side and the potential energy gain on the other side. Indeed the 
condensation energy found as a function of the density of pairs has a maximum at half- 
filling. For example in the weak-coupling BCS regime Econd is proportional to P^ |49] . 
Therefore the condensation energy gain will increase by choosing the number of particles 
in the paired channel as close as possible to half-filling. On the other hand, for a fixed 
total density n, this would reduce or increase the unpaired fermions and consequently 
the potential energy gain, which is maximal for a non-magnetized system since U is 
negative. The competition between these opposite trends eventually determines the 
value of the magnetization in equilibrium, which is finite and rather small at this value 
of the coupling (see inset in Fig. [?]). At half- filling no condensation energy gain can be 
achieved by creating a density imbalance between the superfluid pairs and the unpaired 
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Figure 8. (Color online) c-SF order parameter P (green triangles), magnetization m 
(red squares), average triple occupancy t = {nin2n^) (orange circles) and difference 
between double occupancies in different channels di2 — rfia (blue diamonds) in the c-SF 
phase, plotted as a function of the interaction |J7/VF| for n = 1 and T = {) for the cubic 
lattice in Z) = 3. (Unconstrained, i.e. F = 0) 



fermions since the condensation energy is already maximal. Therefore the spontaneous 
symmetry breaking in the color superfluid phase does not result necessarily in a density 
imbalance, which is however triggered by a condensation energy gain for every density 
deviation from the particle- hole symmetric point. 

We now consider the same system for fixed total density n = 1 and study the ground 
state properties as a function of the interaction strength \U\ (see Figjs]). For weak 
interactions the system is in a c-SF phase. Upon increasing \U\, the order parameter 
P first increases and then shows the dome shape at intermediate couplings which we 
already observed for the half-filled case. Away from the half-filling, the value where P 
reaches its maximum is shifted to lower values of the interaction strength. The triple 
occupancy t, on the other hand monotonically increases with \U\. Interestingly the 
magnetization m{U) has a non-monotonic behavior. At weak-coupling, magnetization 
m{U) grows with increase of the interaction strength. For increasing coupling, m 
has a maximum and then decreases for larger \U\, indicating a non-trivial evolution 
due to competition between the condensation energy and the potential energies for 
increasing attraction. The spontaneous breaking of the SU (3)-symmetry is also well 
visible in the behavior of the double occupancies. Indeed in the c-SF for n < 1.5 
we find di2 > di^ = ^23. The difference di2 — (^23 is however non-monotonic in the 
coupling and seems to vanish at ^ 0.35. Our interpretation is that beyond this 

point the SU (3)-symmetry is restored and the system undergoes a transition to a Fermi 
liquid trionic phase. Indeed for > 0.35 we did not find any converged solution 

within our DMFT approach, neither for a homogeneous nor for a staggered criterion 
of convergence. This result is compatible with the presence of a macroscopically large 
number of degenerate trionic configurations away from the half-filling. A finite kinetic 
energy for the trions would remove this degeneracy, leading to a trionic Fermi liquid 
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Figure 9. (Color online) Superfluid order parameter on the 2D square lattice for 
different total filling as a function of the interaction strength. We neglect spontaneous 
magnetization in the system. (Unconstrained, i.e. F = 0) 
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Figure 10. (Color online) The quasi-particle weight Z averaged over the Fermi surface 
as a function of the interaction strength \U\. (Unconstrained, i.e. y = 0) 



ground state. This contribution is however beyond the DMFT description of the trionic 
phase where trions are immobile objects. We can address the existence of a Fermi hquid 
trionic phase at strong-couphng using the VMC approach in which we will discuss 
in the following. 

As already mentioned in Sec. [3} we use different trial wavefunctions to study the 
behavior of the system in the weak- (|f/| < 1^/2) and the strong-coupling (|f/| > W j'l) 
regimes. At weak-coupling the magnetization is expected to be very small and we 
can consider the results for the unpolarized system with n\ = = to be a good 
approximation of the real system which is in general polarized. We found indeed that for 
IL'^I < W/2 the system is in the c-SF phase with a finite order parameter P. As shown in 
Fig. [9| we obtain that P{U) has a similar dome shape as in the 3D case. Unfortunately, 
we cannot directly address the trionic transition within this approach since it is expected 
to take place at intermediate coupling where both ansatz wave functions are inaccurate. 
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Figure 11. (Color online) Number of particles in the paired channels ni2 = ni = n2 
(blue circles) and the unpaired channel (red squares) and superfluid order parameter 
P as a function of the 3-body repulsion V for |C/|/W^ = 0.312 and total density n = 0.48 
for the cubic lattice in Z) = 3 at zero temperature. Dashed lines correspond to the 
asymptotic values. 



We can however consider the system in the strong-couphng hmit by using the effective 
trionic Hamihonian of Eq. 28 In this way we can study the Fermi hquid trionic phase 
which we characterize by evaluating the quasiparticle weight, averaged over the 
Fermi surface 

Z = %^^. (44) 

Here is extracted from the jump in the momentum distribution at the Fermi surface, 
which we approximate as 

Zk = n^-^ {nk+Ak^ + nk+Aky) , (45) 

where Ak^ i'^^y) is the translational vector along the x (y) direction in the reciprocal 
lattice. In Fig. 10 we plot Z as a function of interaction strength |?7|/iy. 

By combining DMFT and VMC results we therefore have strong evidence of the 
system undergoing a phase transition from a magnetized color-superfluid to a trionic 
Fermi liquid phase at strong-coupling, when the density is far enough from the particle- 
hole symmetric point. 



5. Results: Constrained System {V = oo) 

As referred to in the introduction, actual laboratory implementations of the model under 
investigation using ultracold gases are often affected with three-body losses, which are 
not Pauli suppressed as in the SU{2) case. As discussed in Ref. [H], the three-body loss 
rate 73 shows a strong dependence on the applied magnetic field. Therefore the results 
presented in the previous section essentially apply to the case of cold gases only whenever 
three-body losses are negligible, i.e. 73 <C J, U. In the general case, in order to model 
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Figure 12. (Color online) Number of particles for the paired channels ni2 = ni = 712 
(blue/dark circles) and the unpaired channel (red squares), c-SF order parameter P 
(green triangles) and total double occupancy d = di2 + rfis + ^23 (violet/light circles) 
calculated within DMFT as a function of the interaction strength |f7|/VF for T = 0, 
V ~ 8QW and n = 0.48 (cubic lattice). The dashed green line corresponds to the 
asymptotic value of the superfluid order parameter in the atomic limit Poo, while 
the dotted blue line corresponds to the asymptotic value of the particle density in 
the paired channel, which is also equal to the asymptotic value of the total double 
occupancy. (Constrained case, V ~ SOVl^) 




the system in presence of three-body losses, one needs a non-equihbrium formulation 
where the number of particles is not conserved. However, as shown in Ref. |10], in 
the regime of strong losses 73 ^ J, U, the probability of having triply occupied sites 
vanishes and the system can still be described using a Hamiltonian formulation with 
a dynamically-generated three-body constraint. To take it into account in our DMFT 
formalism, we introduce a three-body repulsion with V = 00. Within VMC we directly 
project triply occupied sites out of the Hilbert space. We stress that finite values of 
V do not correspond to real systems with moderately large 73 since then real losses 
occur and a purely Hamiltonian description does not apply any more; only the limits 
73 ^ J, f/ and 73 ^ J, f/ lend themselves to an effective Hamiltonian formulation. 



5.1. Ground State Properties 

In order to address how the system approaches the constrained regime with increasing 
V, we first used DMFT to study the ground-state properties of the model in 3D as a 
function of the three-body interaction V for a fixed value of the total density n = 0.48 
and the two-body attraction 1^71 /VF = 0.3125. We found that the average number 
of triply occupied sites t = {nin2n^) (not shown) vanishes very fast with increasing 
V. The SF order parameter P and the densities in the paired and unpaired channels 



approach their asymptotic values already for V ~ SVTor V ~ 10|f/|, as shown in Fig. 11 
Therefore, we assume that we can safely consider the system to be in the constrained 
regime whenever V is chosen to be much larger than this value. 
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Figure 13. Effect of the magnetization for total density n — 0.48 on the 2D square 
lattice with 50 lattice sites. In particular we plot AE{ni) — E{m) — E{Q) as a function 
of magnetization m = ni2 — (ni2 = ni — 71,2) for different values of the interaction 
strength U . Calculations are performed using the VMC method with a strong coupling 
ansatz. In the inset we plot /S.E as a function of the interaction strength \U\ for the 
fully magnetized c-SF phase. (Constrained case) 



Both the densities rig- and the superfluid order parameter P are strongly affected 
by the three-body interaction (see Fig. 11). For this value of the interaction, P and 
m are strongly suppressed by the three-body repulsion, even though both eventually 
saturate to a finite value for large enough V . However, as shown below, this suppression 
of the magnetization and SF properties is specific to the weak-coupling regime and for 
larger values of |?7| both the SF order parameter P and the magnetization m are instead 
strongly enhanced in the presence of large V . 

We now investigate the constrained case (setting V = lOOOJ ~ SOW within the 
DMFT approach) where the total density is fixed as above to n = 0.48. Large values 
of the density imply an increase of the probability of real losses over a finite interval of 
time. Therefore we restrict ourselves to a relatively low density which is meant to be 
representative of a possible experimental setup. 

We study the evolution of the ground state of the system in 2D and 3D as a function 
of the two-body interaction strength U. DMFT results in Fig. 12 show that in the three- 
dimensional system the trionic phase at strong coupling is completely suppressed by the 
three-body constraint and the ground state is found to be always a color superfluid 
for any value of the attraction. This remaining c-SF phase shows however a very 
peculiar behavior of the magnetization m as a function of the attraction U. Indeed 
the magnetization m = ni2 — (?^i2 = fii = 712) steadily increases for increasing 
interaction and ^ {m ^ ni2 ~ n/2) already for U ~ 12 J = W. 

Our explanation is that the three-body constraint strongly affects the energetic 
balance within the c-SF phase. Indeed, in the absence of V the magnetization was 
shown to be non-monotonic and to vanish in the S'?7(3)-symmetric trionic phase at 
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strong-coupling. Now instead in the same limit the fully polarized c-SF system has a 
smaller ground state energy for fixed total density n. This result is fully confirmed by 
the VMC data for the 2D square lattice. As shown in the next section, combining these 
results essentially implies that a globally homogeneous phase with m = is unstable in 
the thermodynamic limit with respect to domain formation whenever the global particle 
number in each species N^r = Yli ^i,o- is conserved. By using the canonical ensemble 
approach of VMC, we can indeed address also metastable phases and study the effect 
on the energy of a finite magnetization for fixed total density n = 0.48. In particular 
we study the energy difference between the magnetized system and the unpolarized one 
with the same n, i.e. AE{m) = E{m) — E{0). Results shown in Fig. 13 indicate that 



at strong-coupling the energy decreases for increasing magnetization and the minimum 
in the ground state energy corresponds to the fully polarized system. In the inset of 



Fig. 13, we show AE as a function of the interaction strength for the fully polarized c-SF 
at strong-coupling, which decreases as AE ~ l/|t/|. We also investigated the system 
in the weak-coupling regime, where our calculation shows that AE{m) has a minimum 
for very small values of the magnetization (not shown). This indicates that also in 2D 
the c-SF ground state at weak-coupling is partially magnetized, in complete agreement 
with the three-dimensional results. 
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Figure 14. (a) The supcrfluid order parameter and (b) condensation energy for 2D 
square lattice for different total fillings as a function of the interaction strength. For 
weak coupling we approximate that the system is not magnetized (green lines and 
circles), while for strong coupling we assume that the system is fully polarized, i.e. 
contains only pairs (dashed blue line and squares) . The dotted line corresponds to the 
superfluid order parameter in the atomic limit Poo (Constrained case) 



Within DMFT the order parameter P in the c-SF ground state shown in Fig. 12 
is also increasing with \U\ and saturates at strong coupling to a finite value, which we 
found to be in agreement with the asymptotic value in the atomic limit for the SU (2) 
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symmetric case [38] 

Poo= lim P(U) = 1^/71(2 -n). (46) 

The total number of double occupancies d is also an increasing function of |f/| and 
saturates for very large \U\ to the value ni2 = n/2 as in the strong coupling limit for 
the SU (2) symmetric system. This means that in the ground state the strong coupling 
limit of the SU{3) model is indistinguishable from the SU{2) case for the same total 
density n and two-body interaction U. As we will show in the next subsection, this is 
not any more true if we consider instead finite temperatures. 

Similar considerations on the superfiuid properties in the ground state apply to 
the two-dimensional case studied within the VMC technique. As the magnetization 
in the weak-coupling regime is very small, we approximated it to zero and consider 
an unpolarized system within the weak-coupling ansatz, while at strong-coupling we 
directly consider the system as fully polarized, i.e. containing only pairs. As visible 



in Fig. 14, P shows a similar behavior to the 3D case. Indeed at weak-coupling both, 
DMFT and VMC, show a BCS exponential behavior in the coupling, while at strong- 
coupling P converges to a constant. 

Within VMC we also studied the condensation energy as explained in Sec. [3j Fig. 
[Mj p shows that the condensation energy first increases with the interaction strength U 
as expected in BCS theory, while it decreases as 1/U at strong-coupling as expected in 
the EEC limit for the SU{2) case [38]. Despite the fact that we cannot reliably address 
the intermediate region, there are also indications that the condensation energy has a 
maximum in this region. 

5.2. Finite temperatures 

We also investigated finite-temperatures properties for the three-dimensional case using 



DMFT. In Fig. 15, we show the evolution at finite temperature T of the SF order 
parameter P and of the magnetization m at fixed values of the interaction U. At low 
temperatures, the system is superfiuid and the magnetization finite. With increase of 
the temperature, both P and m decrease and then vanish simultaneously at the critical 
temperature T = Tc{U). This clearly reflects the close connection between superfiuid 
properties and magnetism in the SU (3)-symmetric case and is markedly different from 
the strongly asymmetric case which we studied in Ref. [12], where the density imbalance 
survives well above the critical temperature. 

It is however remarkable that for \U\ > Um ~ W, m{T) and P{T) clearly show in 



Fig. 15 the existence of a plateau at finite T, indicating that the system stays in practice 
fully polarized in a finite range of temperatures. This allows us to define operatively a 
second temperature Tp{U) below which the system is fully polarized, while for T > Tp 
instead the magnetization decreases and eventually vanishes at Tc. 



We summarize these results in the phase diagram in Fig. 16 Inside the region 



marked in orange (|f/| > Um and T < Tp) the system is fully polarized and therefore 
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Figure 15. (a) Superconducting order parameter and (b) magnetization as a function 
of temperature for different values of the interaction strength U. (Constrained case, 
V ~ SOW) 
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Figure 16. (Color online) Phase diagram of the model on the cubic lattice with three- 
body constraint. The solid blue line separates normal and color superfluid phases. 
Below the dashed orange line the system is fully polarized. The dotted black line 
describes the strong coupling behavior of the critical temperature and is obtained by 
a fitting procedure. 
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identical to the SU (2) superfluid case. As we will see in the next section, in a canonical 
ensemble where the total number of particles A'^,^ of each species is fixed, this analogy 
is not any more correct and we have to invoke the presence of domain formation to 
reconcile these findings with the global number conservation in each channel. Outside 
this region and below Tc (solid blue line in Fig. 16), the c-SF is partially magnetized 
and therefore intrinsically different from the case with only two species. This is also 
visible in the behavior of the critical temperature where the SU (3)-symmetry is restored 
in the normal phase. We found indeed that the critical temperature first increases with 
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the interaction strength \U\, similarly to the SU{2) case. Then for \U\ = Um, the 
critical temperature suddenly changes trend and for larger |f/| a power-law decrease 
Tc oc 1/\U\ occurs as shown in Fig. 16 In the SU{2) symmetric case this power-law 
behavior only appears for very large \U\ (bosonic limit) [38], while in the SU{3) case 
this regime occurs immediately for \U\ > Um- The smooth crossover in T^iU) and 
the maximum of in the critical temperature characteristic of the SU{2) case, here are 
replaced by a cusp at |f/| = Um, which marks the abrupt transition from one regime to 
the other. 



6. Domain Formation 



One of the main results of this work is the close connection between superfluidity and 
magnetization in the c-SF phase. Indeed we found that in the c-SF phase, away from 
the particle- hole symmetric point, the magnetization is always non-zero. On the other 
hand ultracold gas experiments are usually performed under conditions where the global 
number of particles = J2i ^i,o- in each hyperfine state is conserved, provided spin flip 
processes are suppressed. The aim of this section is to show that domain formation 
provides a way to reconcile our findings with these circumstances. In particular, 
combining DMFT and VMC findings, we will show that a globally homogeneous c- 
SF phase is unstable with respect to formation of domains with different c-SF phases in 
the thermodynamic limit. 

To be more specific, we will consider the case when the global numbers of particles 
in each species are the same, i.e. A^i = A^2 = ^3 = N/3, at T = 0, though the 
discussion can be easily generalized to other cases. The simplest solution compatible 
with N„ = N/3 is clearly a non-polarized c-SF phase with energy E^om per lattice 
site. This phase is actually unstable and therefore not accessible in a grand canonical 
approach like DMFT, where we fix the global chemical potential fi and calculate the 
particle densities n„ as an output. Since, as shown in Sec. |4] and Sec. |5| the 
system is spontaneously magnetized in the color superfiuid phase out of half-filling, 
there is no way to reconcile the DMFT result with the global constraint N„ = N/3 
assuming the presence of a single homogeneous phase. The VMC approach, on the 
other hand, operates in the canonical ensemble, and it can be used to estimate the 
ground state energy per lattice site for specific trial configurations. For the homogeneous 
configuration, we have Ehom = E{m = 0)„, where n = N/M and M is the number of 
lattice sites. 

Let us now contrast this situation with the spatially non-uniform scenario in 
which we have many color superfiuid domains in equilibrium. Each of these domains 
corresponds to one of the solutions obtained above, and therefore this phenomenon 
can be seen as a special form of phase separation. For two or more phases to be in 
thermodynamic equilibrium with each other at T = 0, they need to have the same 
value of the grand potential per lattice site Q = E — fm for the same given value of 
the chemical potential n, while the onsite density of particles for each species can be 
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(a) (b) (c) 

Figure 17. (Color online) Schematic picture of the phases of a S'C/(3)-symmetric 
mixture of three-species fermions for the total particle numbers in each species 
Na — N/3 away from half filling, (a) visualizes the ground state configuration at 
weak to intermediate coupling in both the unconstrained {V = 0) and the constrained 
{V — oo) case. Irrespective of the presence of the constraint, a finite magnetization 
points at domain formation in experiments with fixed (see text); a specific example 
of a phase-separated configuration is plotted. Increasing the attraction strength reveals 
substantial differences between the two cases: (b) In the constrained case, domain 
formation persists to strong coupling, in parallel to the 3-component asymmetric 
situation [42j . The unpaired species are expelled from the paired regions, pairing up in 
other spatial domains, (c) In the unconstrained case instead, a spatially homogeneous 
trionic phase emerges [501 lU ■ 



different in the different phases. 

Possible candidate phases for the system considered in this paper are suggested by 
the underlying S't/(3)-symmetry. Indeed if we consider c-SF solutions corresponding 
to different gauge fixing, i.e. with pairing in different channels, they will have the 
same total onsite density n and therefore the same energy and grand potential, since 
they correspond to different realizations of the spontaneously broken symmetry. If we 
consider for simplicity only the three solutions with pairing between the natural species 
sketched in Fig. [T7| then this mixture of phases has globally the same number of particles 
iVo- = N/3 in each hyperfine state whenever we choose the fraction of each phase in the 
mixture to be a = 1/3 and n = N/M in each domain. In fact in each domain we have the 
same densities Hp in the paired channel and for the unpaired fermions, even though 
they involve different species in different domains. This scenario is therefore compatible 
with the global number constraint N^r = N/3 and we can compare its energy with the 
energy Ehom of the globally homogeneous c-SF phase. The VMC calculations reported 
in Fig. [13] clearly indicate that for a fixed onsite density n, the ground state energy per 
lattice site is lower by having a finite magnetization, i.e. E{m)n < E{0)n and therefore 
Ehom > Ephase- separated = « Ya=i = ^i^) ^nd Ei = E{m) is the energy per lattice 
site in the i-th domain. Thus a globally homogeneous c-SF phase has higher energy 
than a mixture of polarized domains with the same A^^- and is therefore unstable with 
respect to phase separation. 



It should be noted however, that the configuration sketched in Fig. 17 only 



Magnetism and domain formation in SU (3) -symmetric multi-species Fermi mixtures?)! 

represents the simplest possible scenario compatible with the global boundary conditions 
= N/3. Indeed in the 5'f/(3)-symmetric case we have continuous set of equivalent 
solutions, since solutions obtained continuously rotating the pairing state from 1-2 to 
a generic linear combination of species have the same energy and are therefore equally 
good candidates for the state with domain formation. Moreover, it is well known that 
having a continuous symmetry breaking is intrinsically different from the discrete case, 
because of the presence of Goldstone modes [T7j. In large but finite systems, the surface 
energy at the interface between domains, which is negligible in the thermodynamic limit, 
will become relevant. On one hand a continuous symmetry breaking allows the system to 
reduce the surface energy cost through an arbitrarily small change of the order parameter 
from domain to domain, pointing toward a scenario where a large number of domains 
is preferable in real systems. On the other hand, when the system is finite, increasing 
the number of domains decreases their extension, reducing the bulk contribution which 
eventually defines number and size of the domains at equilibrium. Based on our current 
approaches, we cannot address the issue of what is the real domain configuration in a 
finite system, neither the question if different scenarios with microscopical modulations 
of the SF order parameter take place Similar conclusions concerning the 

emergence of domain formation in the c-SF phase have been already drawn in [20 11211137] 
and also in a very recent work [IH], which addresses the same system in continuum space. 

In real experiments both finite-size effects and inhomogeneities due to the trapping 
potential could play an important role in the actual realization of the presented scenario. 
Furthermore, as the S'?7(3)-symmetry in the cold atomic systems is not fundamental but 
arises as a consequence of fine-tuning of the interaction parameters, imperfections will 
also arise from slight asymmetries in these parameters. We have shown before [12] that 
in the strongly asymmetric limit, phase separation is a very robust phenomenon. We 
may therefore conjecture that interaction parameter asymmetries favor this scenario. 

The combination of the findings in the present paper on the SU{3) case with 
those on the strongly-asymmetric case in [12] suggests that phase-separation in globally 
balanced mixtures is a quite general feature of three-species Fermi mixtures. However, 
the phases involved are in general different in different setups. In the strongly- 
asymmetric case in presence of a three-body constraint, the color superfluid phase 
undergo a spatial separation in superfluid dimers and unpaired fermions [12]. In this 
case, the presence of the constraint is crucial to the phase-separation phenomenon, as 
testified by its survival well above the critical temperature for the disappearance of the 
superfluid phase [12]. In the fully SU{3) symmetric case instead, the presence of the 
constraint only modifies the nature of the underlying color superfluid phase favoring 
fully polarized domains at strong coupling. The formation of many equivalent color 
superfluid domains can be seen as a special case of phase separation reflecting the SU (3) 
symmetry. In this case the phase separation phenomenon is strongly connected to the 
superfluid and magnetic properties of the color superfluid phase and it is expected to 
disappear at the critical temperature Tc and for the peculiar particle-hole symmetric 
point at half- filling in the unconstrained case. 
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7. Conclusions 

We have studied a SU(3) attractively interacting mixture of three-species fermions in 
a lattice with and without a three-body constraint using dynamical mean-field theory 
{D > 3) and variational Monte Carlo techniques {D = 2). We have investigated both 
ground state properties of the system and the effect of finite temperature and find a 
rich phase diagram. 

For the unconstrained system, we found a phase transition from a color superfiuid 
state to a trionic phase, which shows additional charge density modulation at half-filling. 
The superfiuid order as well as CDW disappear with increasing temperature. 

In the presence of the three-body constraint, the ground state is always superfiuid, 
but for strong interactions \U\ > Um the system becomes fully polarized for fixed total 
density n. It is remarkable that according to our calculations the system stays fully 
polarized in a range of low temperatures. For high temperatures a transition to the 
non-superfiuid SU (3) Fermi liquid phase is found. The critical temperature has a cusp 
precisely at Um- This is in contrast to the S'f/(2)-symmetric case, where a smooth 
crossover in the critical temperature takes place. 

The c-SF phase shows an interesting interplay between superfiuid and magnetic 
properties. Except in the special case of half-filling, the c-SF phase always implies a 
spontaneous magnetization which leads to domain formation in balanced 3-component 
mixture. 
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Appendix A. Derivation of the strong coupling Hamiltonians 

Appendix A.l. Constrained case 

In order to derive a perturbative strong-coupling Hamiltonian for the constrained case 
we make use of the Wolff-Schrieffer transformation [17] 

Hpert = Vne'^He-'^VD (A.l) 

and keep terms up to the second order in J/Uaa'- In the expression above, Vd is the 
projection operator to the Hilbert subspace with fixed numbers of double occupancies 
in each channel {Np, Aj^,A]'^), and e**^ is a unitary transformation defined below. 
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The kinetic energy operator can be split in several contributions, where the subscripts 
indicate the change in the total number of double occupancies (A^d,o = + + 



13N 
d )i 



I.e. 



/Co 



/C-1 



(A.2) 
(A.3) 
(A.4) 



Here = 4^Cj /li,^ = 1 - n^^^ and a ^ a a a. 

We note that whereas JCq preserves the total double occupancy N^fl, it contains two 
different types of terms: (i) terms that also preserve double occupancy in each channel 
N^'^' ( /Cg part) and (ii) terms that change the double occupancy in two different channels 
such that the total double occupancy stays unchanged (/Cg part). Thus, we can write 

/Co = /Co" + /C^. (A.5) 

We can also decompose the operators that change the total number of double 
occupancies into 

/Ci = /C}2 + /Cf + /CJ^ (A.6) 
/C_i = ]Cl\ + JC^\ + ]Ci\ , (A.7) 

where the superscripts give the type of double occupancies that are being created or 
destroyed. The canonical transformation can be written as an expansion to the second 
order 

(A.8) 



where H = /Cq + /Ci + /C_i + V and V = Yl U^^'Tii^^ni ^ji . Then, we choose 

i,(T<(7' 

= E {^(^r' - ^-'O + ([/cr',/Co] + [/c^r,/Co 



.(A.9) 



Inserting Eq. (A.9) into the Eq. (A.8) we obtain 



"Hpert — V + /Cq 



(A.IO) 



<T<cr' (t"<(t"' 



2Ua-a'Ucr"a' 



-V 



D 



{icf - /cr/), [/Cf '^"^ V] - [/C^7", V] Vd + o{ 



U2' 



Using the relation [V,/C5.i'] = ±Ufj^„iK'^' and applying the projection Vd, we arrive at 



n 



pert 



a<<7' 



'IP 



(A.ll) 
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Notice that most of the terms in the commutator become zero leaving only the correlated 
hopping terms. 



In order to write Eq. ( A.l 1 ) in a more practical way, we can define double occupancy 
operators as d\^^, = cl^^Uiyhiy and single occupancy operators as f^^ = hiyhiy/c\^^ 
with a ^ a' ^ a" ^ a. With this notation, the perturbative Hamiltonian becomes 



(«jV (j,t);(i,j);a<a' 



J' E T^4,.'J^yflA,<r.' + V + . (A.12) 



For the case where the S'f/(3)-symmetry is restored {Uaa' = U), the perturbative 
Hamiltonian can be written in a compact notation 

H,ert = ^-jJ2 [fih'^ + 4, A,'] - ^ E 4,.AJlA'^ (A- 13) 

J2 J2 j3 

where the double occupancy operator is now defined as = c|g.(/ii,o-''^i,o-" + hiyiUiy). 
Appendix A. 2. Unconstrained case 

Without the 3-body constraint three fermions with different hyperfine states can occupy 
the same lattice site and we expect them to form trionic bound states at sufficiently 
strong coupling. 

According to perturbation theory up to third order we could have two different 
contributions: (i) one of the fermions hops to one of the neighboring sites and returns 
back to the original site (second order perturbation), (ii) all three fermions hop to the 
same nearest neighbor site (third order perturbation). As we show below, due to the 
first process there is an effective interaction between trions on nearest neighbor sites. 
Also due to this process the onsite energy has to be renormalized. The second process 
(ii) describes the hopping of a local trion to a neighboring site. 

The energy gain due to virtual processes, when one of the fermions is hopping to 
a nearest neighboring site and returning back, can be easily determined within second- 
order perturbation theory 



Eto — Ei, 

l.(T " 



where Yl'i denotes summation only over the nearest neighbors of the trion. Here |to) 
describes a local trionic state at lattice site 0, while by \ia) we define a state where site 
i is occupied by a fermion with spin a, while two other fermions stay in the lattice site 
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0. One can easily calculate that |(z(j|?^|to)P = and Et^^ — Eia- = U^a' + Uaa", where 
a ^ a' a" a. So we obtain 

where z is the number of the nearest neighbor lattice sites. 

The calculation above assumes that neighboring sites of a trion are not occupied. 
If one of the neighboring sites is occupied by another trion, then the energy gain per 
trion is given by 

The effective interaction between two trions on neighboring sites is therefore 

J2 J2 J2 



For the SU (3)-symmetric case this expression is simplified and we obtain 

Therefore the nearest neighbor interaction between trions is repulsive in the SU{3)- 
symmetric case. 

The next step is to calculate the effective hopping of the trions. For this purpose 
one has to use third order perturbation theory 

{to\n\a){a\n\aa'){aa'\n\t^) 

{E,-E^)[E,-E^^,) ■ ^^-''^ 

Here \to) and \ti) define local trions on lattice site and the neighboring lattice site 1 
respectively, \a) defines a state where a fermion with spin a occupies the lattice site 1, 
and two other fermions are occupying the lattice site 0. Conversely Icrcr') defines a state 
where two fermions with spins a and a' occupy the lattice site 1. On the lattice site 
we have only a fermion with spin a" 7^ a, a'. For any a and a' the matrix elements 
are given by (to|'H|cr) = {alT-Llaa') = (aa'lT-Llti) = —J, Et^^ — E„ = Uaa' + Uaa" and 
Et^ — Ea-a' = Ua^ii + Ua^'a", whcre 0", 0"' and a" are three different hyperfine-spins. 
So we obtain 

where a, a' and a" are different from each other in the sum. 

In the 5'f/(3)-symmetric case, the expression again simplifies to 

Jeff = (A.21) 



So we obtain the following effective Hamiltonian 

^e// = -Jeff Yl tlt^ + Veff ujuj . (A.22) 

Here t] is the creation operator of a local trion at lattice site i and nj = t\ti is the 
trionic number operator. 
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